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Hard-sphere mixtures provide one a solvable reference system that can be used to improve the density func- 
tional theory of realistic molecular fluids. We show how the Kierlik-Rosinberg's scalar version of the funda- 
mental measure density functional theory of hard spheres [Phys. Rev. A, 42, 3382 (1990)], which presents 
computational advantages with respect to the original Rosenfeld's vectorial formulation or its extensions, can 
be implemented and minimized in three dimensions to describe fluid mixtures in complex environments. This 
implementation is used as a basis for defining a molecular density functional theory of water around molecular 
hydrophobic solutes of arbitrary shape. 



I. INTRODUCTION 

The numerical methods that have emerged in the sec- 
ond part of the last century from liquid-state theories^l, 
including integral equation theory in the interaction- 
site^Ml or moleculaJ^"^ picture, classical density func- 
tional theory (DFTjpHHJ or classical fields theorjffaHUl 
have become methods of choice for many physical chem- 
istry or chemical engineering applications^^!]. They 
can yield reliable predictions for both the microscopic 
structure and the thermodynamic properties of molecu- 
lar fluids in bulk, intcrfacial, or confined conditions at a 
much more modest computational cost than molecular- 
dynamics or Monte-Carlo simulations. A current chal- 
lenge concerns their implementation in three dimensions 
in order to describe molecular liquids, solutions, and 
mixtures in complex environments such as atomistically- 
resolved solid interfaces or biomolecular media. There 
have been a number of recent efforts in that direc- 
tion using 3D-RIS1VP^ ^, m olecular density functional 
theorySHlEI] lattice fielcP^or Gaussian fiekP theories. 

For condensed homogeneous and inhomogeneous flu- 
ids, the hard-sphere (HS) model plays a central role. It 
provides not only a good physical representation of col- 
loidal dispersions where the range of inter-particle at- 
traction is typically much smaller than the particle size, 
but also an invaluable reference system for studying the 
properties of simple liquids where the structure is pre- 
dominantly determined by the short-ranged repulsion. 
In this respect, following t he pre cusor work of Percus 
for one dimensional systems^! ^ and weighted density 
ideas by Tarazona and EvanstiMiiJ the Rosenfeld's deriva- 
tion in 1989 of a quasi-exact DFT for inhomogeneous 
hard-sphere mixtures in three dimensions, the fundamen- 
tal measure theory (FMT), constitutes a major advent 
of modern statistical mechanics.^ Several extensions or 
variants of Rosenfeld's FMT were proposed subsequently. 
One year later, an alternative, scalar rather than vec- 
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torial formulati on of FMT was derived by Kierlik and 
Rosinbcrg (KR^Mill this formulation (which was later 
shown to be mathematically equivalent to the Rosenfeld's 
original functional^) will be the focus of the present 
work. On the other hand, it was soon realized that FMT 
in its original form was able to describe very accurately 
fluids at interfaces but showed serious limitations for 
the description of the solid phase and liquid-solid tran- 
sition, or that of highly confined fluids. Several success- 
ful solutions were proposed in the following two decades 
to extend FMT to the solioP2HlZl ; including a tensorial 
correction to the vectorial formulation that is able to 
satisfy various dimensional crossover requirements^^ 
. Besides, an extension of the Rosenfeld's vectorial ver- 
sion, based on the Mansoori-Carnahan-Starling-Leland 
(MCSL) hard-sphere equation of state rather than the 
Percus- Yevick (PY) one, was proposed independently by 
Roth et alPO (the so-called White-Bear (WB) version) 
and Wu et alP^ (modified fundamental measure theory, 
MFMT). The WB version was made compatible with 
the Tarazona's tensorial extensions to describe crystalline 
phases^. For recent reviews of FMT and of the twenty 
years of efforts that have followe d to im prove on the ini- 
tial Rosenfeld's proposal, see Ref d 60 * 61 l 

The existence of a hard-sphere reference for the DFT of 
classical fluids has promoted recently a great deal of ap- 
plications of thi s app roach to the study of atomic- like and 
polymeric fluids 62 63]i , as well as simplified models of aque- 
ous solutions^ or ionic liquids^ 6 ^, both at interfaces or 
in confinement. Most of those studies arc limited how- 
ever to planar, cylindrical or spherical symmetries. To 
date, there have been few applications of classical DFT 
in three dimensions, in orde r to cop e with fluids in com- 
plex molecular environment j 3 ° l 3 H SZI. I n particular only a 
few 3D implementation of FMT have been described in 
the literature^H 6 ^. All of them are based on the origi- 
nal Rosenfeld's vectorial formulation. This formulation 
involves four scalar weigthed densities and two vecto- 
rial ones, making a total of ten weighted density compo- 
nents to be handled. In this paper, we propose the first 
3D implementation of the Kierlik and Rosinberg scalar 
version of FMT (KR-FMT) , using the Percus- Yevick or 
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Carnahan-Starling variants that are both described in the 
original KR paper. This formulation involves only four 
scalar weighted densities which, as will be seen, leads to a 
substantial computer speedup for three-dimensional ap- 
plications with respect to the vectorial version, especially 
for multicomponent systems. 

The needs for a three-dimensional implementations of 
FMT are broad and go beyond the hard-sphere model per 
sc. Obviously, the hard-sphere functional can be used as 
a reference for representing the hard core interactions 
in a molecular system, and other interactions such as 
Lennard- Jones or Yukawa attraction and Coulombic in- 
teractions can be added in the functional as a mean-field 
perturbation in order to model realistic molecular fluids, 
including water and ionic solutions^. Such functionals 
deserve further developments and applications to the case 
of complex molecular environments. Our goal is here an- 
other. We have constructed recently a three-dimensional 
molecular density functional theory (MDFT) approach to 
solvation in molecular liquids that is based on the knowl- 
edge of the angular-dependent direct correlation function 
of the pure solvent (the so-called homogeneous reference 
fluid approximation, HRF^J). This approach amounts 
to a second-order Taylor expansion of the free-energy 
with respect to the density around the homogeneous 
fluid reference. It is connected to the hypernetted chain 
approximation (HNC) in integral equation approaches 
and amounts to incorporate in the functional only the 
homogeneous-solvent two-body correlations. Such ap- 
proximation proved to be accurate for polar molecular 
fluids such as acetonitril^22H25l ; brrt appeared clearly in- 
sufficient in the case of water. In RefP^, we proposed 
to introduce empirical three-body correction terms in- 
spired by the Molinero's water modeP^ that re-introduce 
some missing tetrahedral symmetry in the HRF func- 
tional. Another possible many-body corrections are those 
induced by the hard-core interactions and that can be 
inferred from the hard-sphere FMT functional. Such ap- 
proximation is at the heart of the reference HNC approx- 
imation (RHNC) in integral equation theories-^ and was 
pushed by Rosenfeld for DFT^l ft was invoked recently 
by Zhao et al in a post-treatment of ionic microscopic sol- 
vation profi les by DFT in order to estimate the solvation 
free-energie d 36 * 37 l This is thus the type of corrections 
that we investigate here. We limit ourselves in a first 
step to the solvation of hydrophobic solutes which, be- 
sides yielding simpler functional forms (the solvent angu- 
lar dependence may be omitted) , deserves special studies 
since hydrophobic solvation has been recently at the cen- 
ter of many debates^HZH. Let us state from the beginning 
that we will be dealing in this paper with the hydration 
of microscopic solutes-^l and that, at the present stage, 
macroscopic effects such as dewetting are not intended 
to be contained in the functional. 

This paper is organized as follows. Section 2 recalls 
the basic principles of hard-sphere fundamental measure 
theory in the original Rosenfeld's vectorial formulation 
and in the Kierlik-Rosinberg scalar version. Section 3 



discusses the relative practical merits of the two formu- 
lations, describes the implementation of the KR func- 
tional in three dimensions, as well as a few tests of the 
method. In section 4, the FMT 3D-implementation is 
used to add N-body hard-sphere corrections to a density 
functional for water in order to describe the solvation of 
microscopic hydrophobic solutes. 



II. FUNDAMENTAL MEASURE THEORY: SCALAR 
VERSUS VECTORIAL FORMULATION 

We consider a model fluid mixture composed of N s 
species represented by hard spheres of radius Ri and bulk 
density pf. The fluid is subjected to an external pertur- 
bation, for example a solid interface or a molecular solute 
of arbitrary shape embedded in the fluid, that creates 
for each species i a position-dependent external poten- 
tial Vi(r) and thus an inhomogeneous density Pi{v). The 
grand potential of the perturbed system can be expressed 
as a functional of the inhomogeneous densities and can 
be evaluated relatively to the homogeneous fluid 



n[{ Pi (r)}}=F[{ Pi (r)}}+n[{pn] 



(1) 



Following th e gen eral scheme of classical density func- 
tional theory™^, the functional T({p l {r)}) can be de- 
composed into an ideal, an external and an excess con- 
tribution, according to 

%W}] = ^[{ftW}]+^[{ftW}] 

+ Fexc[{Pi{*)} ~ -^exc [{/?"}] 



dr( Pl {v)-p*) 



(2) 



where 



F id [{ P i{r)}] =k B Tj2 [ dr Pi (r)]xi 

-Pi(r)+P°i 
JW{*(r)}]=£ / drV t (r) Pl (r) 



Pi(r) 



(3) 
(4) 



with ks is the Boltzmann constant and T is the tem- 
perature. jF exc ({pi(r)}) is the excess functional for the 
hard-sphere fluid and p l exc is the bulk excess chemical 
potential of each species defined by 

^ exc ~ — Jp~(7) — \{p>w}={p°} w 

In the fundamental measure theory introduced by 
Rosenfeld'™, the excess functional for the hard-sphere 
fluid can be written in terms of a set of N w weighted 
densities, n a (r): 



•W{Mr)}] = k B T / dr$({n a (r)}) 



(6) 
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with 



by 



n a (r) = J2 J dr' Pi (r')<(r-r') 



Pi(r)*^(r) 



( 7 ) 

where wj,(r) are geometrical weight functions to be de- 
fined below and * indicates the convolution of the mi- 
croscopic densities by those weight functions. The func- 
tional derivative of this excess free energy with respect 
to the densities is given by: 



"' k B Tj2 



dr' 



= k R T 



d$ dn a (r') 
dn a (r') dpi(r) 



= k R T 



„ dn a (r) 



dn a (r' 

( r )> 



(8) 



which appears as the convolution of the partial deriva- 
tives of $ by the weight functions. The equilibrium in- 
homogeneous densities in the presence of the external 
potential Vi (r) are obtained by minimization of the func- 
tional defined above, which is equivalent to solving the 
following Euler-Lagrange equation for all of the species 



— — - = fc B Tln — q +Vi{r + 
fair) \ Pi J $P*(r) 



"Me 



(9) 



In the original Rosenfeld's derivation there are four 
scalar weight functions, uj l a (r), a = 0,1, 2, 3, and two vec- 
torial ones t3i(r), ^(r) per species i that are defined by 



a4(r) = Q(Ri - 



(10) 



wS(r)=47ri2 i wl(r)=47ri2?w5(r) = «(i2 i -r) (11) 



t^(r)=47ri2i£^(r) = - <5(J2i - r) 



(12) 



O(r) denotes the Heaviside function and 8(r) the Dirac 
distribution. The excess free-energy density $ derived by 
Rosenfeld for Eq. [6] is a function of the three position- 
dependent weighted densities, n a (r),a = 0,1,2,3, and 
of the two vectorial ones, ni(r), n 2 (r), which generates 
in the homogeneous limit the Percus-Yevick equation of 
state for hard-sphere mixtures. Starting from the gener- 
alization of the Carnahan-Starling (CS) equation of state 
to mixtures (namely the Mansoori-Carnahan-Starling- 
Leland equation (MCSL)) instead of PY, Roth et aP 1 
and Wu et af^l were later able to obtain a modified ex- 
pression based on the same definition of the weighted 
densities (either called white-bear (WB) version or mod- 
ified FMT version (MFMT)). This modified version of 
FMT takes advantage of the fact that the CS expression 
provides one a better equation of state that PY. 

Ten years before those latest developments, Kierlik and 
Rosinberg were able to derive an alternative version of 
FMT which involv es on ly four scalar weight functions 
u4(r), a = 0, 1, 2, 3PE3 The last two weights are iden- 
tical to Eq. TofTT whereas the first two ones are given 



oji(r) = —S'(R i -r) 



(13) 



W j(r) = ^ S"{Ri -r) + ^ 8'(Ri r) (14) 

Those weight functions appear naturally in the derivation 
as the inverse Fourier transforms of 

4-7T 

uUk) = — (sin(fci?j) - kR l cos(kR i j) 
u\{k) — sin(kRi) 

u Uk) = — (sm(kRi) + kRiCostkRi)) (15) 
2k 

kR. 

Wo(fc) = cos(kRi) H — — sin(kRi) 

Although the main part of the papers by Kierlik and 
Rosinberg relies on a PY expression for the excess free 
energy density 



4P v K] = _ noh(1 _„ 3)+ i^ + _l_ «. , (16) 



the authors do mention in their conclusion that a CS 
(more precisely MCSL) expression could be used instead 



1 n 2 \ i i 

36^^- no I ,n(l "" ;i 
nin 2 



1 — n 3 36n (1 — n 3 ) 2 n 3 



(17) 



They point out the fact that this expression is more pre- 
cise than the PY one, but using it while keeping the ex- 
pression of the weights unchanged leads to thermody- 
namic inconsistencies; those inconsistencies are indeed 
present in the WB or MFMT formulations too. There 
is clearly a trade off to be made between precision and 
theoretical consistency. It was later shown by Phan et 
al. that the Kierlik and Rosinberg's approach is math- 
ematically equivalent to the original vectorial version. 80 
On a practical point of view, however, and especially in 
the perspective of 3D applications, the KR formulation 
is advantageous with respect to the Rosenfeld's formula- 
tion since the number of weighted densities is reduced. 
So is the number of convolutions and thus the number of 
3D-Fast Fourier Transforms (3D-FFT) to be performed. 
This technical point is discussed below in more details. 
Before proceeding, we note again that the functionals 
considered above are well suited to describe inhomoge- 
neous liquids at interfaces or in loose confinement, but 
they are known to fail for crystalline phases or highly 
confined conditions. In thos e cases, various extensions of 
3D-FMT have been devisecPHIElEU. They lie outside 
the scope of the present study. 
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III. IMPLEMENTATION OF THE KIERLIK-ROSINBERG 
FUNCTIONAL IN 3D 

As mentioned in the introduction, there have been a 
few 3D-implementations of FMT proposed in the lit- 
erature, the first one by Frink and Salinger^. All of 
them are based on the Rosenfeld's vectorial formulation. 
Whatever the formulation chosen, however, a natural way 
to solve the FMT equations in 3D is to discretize them 
on a 3D orthorombic grid and to handle the convolutions 
through 3D-FFT's. A typical minimization algorithm re- 
quires one to provide at each minimization step the value 
of the functional and of the functional derivatives for a 
given set of densities {pi(r)}. If we denote by N w the 
number of weight functions to be considered, the numer- 
ical procedure involves 1) to transform the densities in 
Fourier space to {/?i(k)}, 2) to multiply those densities by 
the N w x N s weight functions and sum the products over 
the different species to get the weighted densities (Eq. [7]) , 
3) to transform back the N w weighted densities to real 
space to compute the excess free energies (eq. [6]) and the 
partial derivatives with respect to those weighted densi- 
ties J^-(r), 4) to transform those quantities to k-space 
and, for each species, to multiply them by the weight 
functions and sum, and finally 5) to back transform the 

SF exc 

results to real space to get - e . xc . for all of the species 

5pi(r) 

(Eq.§. 

The whole procedure sums up to 2(N S + N w ) FFT's 
to be performed. N w = 4 for the KR scheme whereas 
N w = 10 in the vectorial formulation (4 scalar weights 
+ (2x3) vectorial components in 3D). For a one com- 
ponent system, one can take advantage in the vectorial 
formulation of the relationships that exist between the 
weights and reduce the number of independent weights 
to be considered to N w = 5. In this case the speedup of 
the scalar versus the vectorial formulation appears rather 
marginal: each cycle requires 10 forward and backward 
FFT's instead of 12, thus a - 20% difference. The re- 
duction of the number of independent weights does not 
apply to multi-component mixtures. The balance thus 
becomes 12 FFT's versus 24 for a two-component system 
and (8 + 2N s ) versus (20 + 2A^) in the general case. The 
expected speedup is thus, in this more generic case, of 
100% and more. 

For those reasons, we propose in this paper the first 
3D-implementation of the Kierlik-Rosinberg's version 
of FMT. We have implemented both the PY and CS 
(MCSL) versions, which only differ in the expression of 
the excess free-energy density, Eq.[l6]or[T7j and the corre- 
sponding partial derivatives with respect to the weighted 
densities. For an arbitrary number of species in the HS 
mixture, the densities are discretized on an orthorombic 
grid of dimension L x x L y x L z . The external potentials 
Vi(r) for every species are first pre-computed and tab- 
ulated. This potential might originate from hard walls, 
or from molecular solutes embedded in the mixture and 
described in terms of site-distributed hard-sphere repul- 



sions or Lennard- Jones interactions. For a given exter- 
nal potential, the FMT functional described by Eq. 2][7 



then minimized using the forward-backward FFT scheme 
described above. The minimization is performed in di- 
rect space with respect to the fictitious "wave-functions" 
^i(r), defined by Pi(r) — ipi{r) 2 , in order to avoid spu- 
rious negative values of the densities that would make 
the logarithm term of the ideal part of the free energy 
functional diverge. As a minimization algorithm, we 
had recourse to the L-BFGS quasi-Newton optimization 
routine 81 which is optimized to handle very large systems 
and requires one, at each step, to provide free-energy 
value and gradients. The gradients are known analyti- 
cally as in eq. [9] 

We first illustrate our FMT implementation for multi- 
component mixtures for the classical test case of a two 
component hard- spher e fluid with radii R\ and i?2 = 3i?i 
near a hard wal l 69 * 72 !, for which reference Monte- Carlo 
calculations are available in the literature^. We show 
in Fig [TJ that the code converges and gives sensible re- 
sults for even a low grid resolution of 3 points per small 
hard-sphere diameter, o\ = 2Ri, in all directions. The 
3D KR-FMT results appear already in excellent agree- 
ment with the simulation data when using a finer res- 
olution of 6 points/<7i in the z-direction while leaving 
the xy-resolution unchanged. A convergence criterium of 
10~ 6 is reached after 10 minimization iterations start- 
ing a a uniform bulk mixture. Such typical conver- 
gence is further illustrated in Fig. [2] for the solvation of 
a neutral benzene molecule, represented by 12 lennard- 
Jones atomic sites in a one-component Lennard- Jones 
fluid modeled by a hard-sphere FMT functional with ra- 
dius R — 1.25 A and at a liquid density po = 0.03328 
particles/A 3 (a simplified representation of water at am- 
bient conditions, see the next section). Starting with 
a guess density p(r) — p exp(— /3V(rj), where V(r) is 
the external potential, it is seen that the convergence is 
basically exponential as a function of the minimization 
step and that typically 10-15 iterations are required to 
get fully converged results. We show in Fig. [2] that the 
required computer time grows linearly with the number 
N g of grid points. Performed on the single processor of a 
standard laptop or desktop computer, a full minimization 
cycle requires between a few seconds for N g — 64 3 and 
a few minutes for N g = 256 3 . For a solvent of the size 
of water (a ~ 3.0A), a typical resolution of 3-4 points/A 
is sufficient to get accurate free-energies and densities. 
With such resolution, one can foresee a possible appli- 
cation of the method to rather large molecular system, 
requiring box sizes up to 100 A. 



IV. APPLICATION TO HYDROPHOBIC SOLVATION 

In continuation to previous works on molecular den- 
sity functional theor y 30 " 3 \ we apply our implementation 
of KR-FMT to improve the MDFT description of molec- 
ular solutes in liquid water at ambient conditions (bulk 
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Distance From Wall [z/2R ] 



FIG. 1. Reduced density profiles 8Rip(z) versus the distance 
from the wall z/2R\ of a binary hard-sphere mixture near 
a hard wall at diameter ratio R2/R1 = 3 and bulk reduced 
densities 0.0260 and 0.0104. The lines represent the 3D-FMT 
results using 3 or 6 grid points per hard-sphere diameter in the 
z-direction (dashed and solid lines, respectively). The black 
dots arc the Monte-Carlo reference simulation data from Tan 
et alP. 
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FIG. 2. Top: Typical plot of the free energy difference be- 
tween two successive steps (normalized by the initial energy) 
versus L-BFGS minimization-step number (Here a benzene 
molecule in a one-component HS reference fluid modeling SPC 
water; see Fig. B. The inlet represents the same with a log- 
arithmic scale in ordinates. Bottom: CPU time per mini- 
mization step versus number of 3D-grid points. The circle 
correspond in increasing order to N= 32, 64, 128, and 256. 



density p = 0.03328 molecules/A 3 ). In MDFT, the sol- 
vent response to the solute external field is described in 
terms of a functional of the inhomogeneous position and 
orientation density, p(r, f2). If the solute is modeled as a 
purely hydrophobic entity, bearing no electrostatic mul- 
tipole, and if the solvent position-orientation coupling is 
neglected (which is a good approximation for water) , the 
angular dependence can be omitted, and the functional 
can be expressed in terms of the isotropic number density, 
p(r) = fcMp(r,n): 



F[p(r)} = k B T J dr 



p(r) In 



Po 



-P(r) 



Po 



drV(r)p(r) 



~HT / dr *'Mr)Ap(r')c s (|r-r'| 
-Mp(r)} 



Po) 
(18) 



with Ap(r) — p(r) — pq. V(r) is the external poten- 
tial to be denned below. The last two terms represent 
the excess free energy, J-" exc [p(r)], which is decomposed 
into a homogeneous reference fluid (HRF) term, involv- 
ing the isotropic direct correlation function of the pure 
solvent at the density po, cg(r; po), and a correction term 
or "bridge" term (in reference to integral equation the- 
ory), which is basically unknown, but can be formally 
expanded in terms of the pure solvent three-body,..., N- 
body direct correlation functions. Th e appro ximation 
that we propose here is rather standar d 36 ^ 72 ! 83 ! and con- 
sists in replacing the unknown bridge for liquid water by 
the exact bridge of an equivalent hard-sphere fluid 



Fb\p(t)] 



^exc [Po] 



fexc 



drAp(r) 



k B T 
2 



J drdr'Ap(r) Ap(r') cf 5 (|r - r'| 



The first three terms represent the one-component hard- 
sphere KR-FMT excess functional defined in the previ- 
ous section and the associated chemical potential yielding 
equilibrium at p(r) = pg. The fourth term involves the 
direct correlation function of the HS fluid at the same 
density, i.c 



MS 1 



\Po) 



5 2 ^[P\ 



(20) 



Sp{r)Sp(r'y p{r)=p °- 
This function can be easily obtained in Fourier space as 
<9 2 $ 



E 



dn a dnp 



(K})w«(*H(*) (21) 



where {n®} represent the weighted densities for a uni- 
form fluid of density po and the U) a> p(k) are the weights 
of eq. 
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The second derivatives have to be taken for 
the PY or CS functions of eqs [16] or [TT] the correspond- 
ing functions are reported in the Appendix. Note that 



defined as in eq. 19 ^^[^(r)] carries an expansion in Ap 



6 



of order 3 and higher that corrects the second order ex- 
pansion of the excess free energy in eq. 
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In this approach, two elements should be further de- 
fined. First the direct correlation of water. In princi- 
ple it can be extracted from the experimental oxygen- 
oxygen structure factor. Since our further comparison 
will be with respect to molecular dynamics simulations 
carried out with the SPC water model, we have com- 
puted cs(r',po) for this model. To do so we have re- 
generated the well-known oxygen-oxygen pair distribu- 
tion function by carrying out molecular dynamics sim- 
ulations with 4096 water molecules and a box size of 
~ 50A; see Fig. [3] The corresponding direct correlation 
function can be deduced by solving the Ornstein-Zernike 
equation. This can be done quite naturally in Fourier- 
space^. To avoid the numerical problems that occur in 
this case at small k-values, we have used instead the di- 
rect space method of BaxteJ^ combined with the vari- 
ational method of Dixon and Hutchinson^ see RefP^ 
for details. This method imposes that the direct corre- 
lation function vanishes beyond a cut-off value that we 
set to R c — 8.7A The corresponding function cs{r;po) 
is plotted in Fig. [3] A second necessary ingredient is to 
fix the radius R of the equivalent hard-sphere fluid. A 
natural value is around R ~ 1.25 A that corresponds to 
the hard core, i.e. the region where hs(r) — in Fig. [3] 
This choice can be confirmed by a striking fact noticed by 
Chandler and Varilly in RefP^ (see their figure 6): the 
statistics of spontaneous empty cavities in SPC water, 
that they tightly link in their Gaussian field analysis to 
the mechanism of hydrophobic solvation, are quite simi- 
lar to those obtained for a hard-sphere liquid at a reduced 
density p* = 8p Q R 3 = 0.5, yielding R ~ 1.25A at the wa- 
ter density. Small variations around that value can be 
conceived and for reasons described below we were led to 
choose R = 1.27A. 

With the previous elements in hand, the functional 
of eqs |18|19| can be minimized in the external Lennard- 
Jones potential field imposed by a molecular solute 
placed at the center of a cubic box. It is defined as 







6" 








\\T-Ti\J 







(22) 



where the 's stand for the positions of the solute atomic 
sites and cr W i,e W i are the site- water Lennard- Jones pa- 
rameters (using Lorentz-Bcrthclot mixing rules and the 
SPC parameters for water). 

For the the KR-FMT excess free-energy terms, we 
have used either the PY or CS version, eqs [16] and 



17 (with very little influence of this choice on the re- 
su ts, as will be seen). In addition to those terms that 
can be handled as described in the previous section, 
one needs to compute the eg and Cg S quadratic terms; 
since they appear as convolutions, this is easily done by 
forward-backward Fourier transform as in a regular HRF 
approximation^!. The procedure is illustrated in Fig. |4j 
representing the three-dimensional density obtained by 




FIG. 3. Left: Oxygen-oxygen isotropic pair distribution func- 
tion of SPC water, hs(r), computed by molecular dynamics. 
Right: corresponding direct correlation function, cs(r), ob- 
tained by Ornstein-Zernike inversion using the Baxter direct- 
space method. 



minimization around a benzene molecule (Lennard- Jones 
parameters from RefP^, no electrostatics). In Fig. [SJ we 
compare the molecule site-water oxygen (C-O w and H- 
O w ) pair distribution functions obtained by DFT with 
and without the bridge term of eq. [19] to the same quan- 
tities generated by molecular dynamics simulations (one 
solute and 512 SPC water molecules). Two features 
should be noted. First the bridge term turns out to have 
little influence on the overall microscopic structure, al- 
though it has a noticeable influence on the computed 
solvation free energies (see below). Secondly, the agree- 
ment to MD can be qualified as quite satisfactory, de- 
spite a disagreement in the shape of the first peak for 
C-O lu . The same type of comparison is drawn in Fig. [6] 
for propane. The geometry and parameters of Ashbaugh 
et alP^I were used with a unified description of the CH 2 
and CH3 groups. Again the agreement with the MD re- 
sults is quite good, with a slight underestimation of the 
first peak width for CH^-O^ . The influence of the bridge 
term on the structure remains marginal. 

Finally Table 1 compares the solvation free energies ob- 
tained for a series of rare-gas atoms and alcanes molecules 
to the MD results reported by Guillot and Guissani us- 
ing a particle insertion method^ or Ashbaugh et al. us- 
ing thermodynamic integration techniques 87 . It is seen 
that the straight application of the HRF approximation, 
J~b[p\ — 0, systematically overestimate the solvation free- 
energies. In Fig. [7] we display the free-energy energy of 
methane obtained when adding the hard-sphere bridge 



term of eq. 19 and varying the hard-sphere radius around 
1.25A. It can be observed that the computed free-energy 
decreases steadily with increasing R (whereas the micro- 
scopic water structure remains unaffected by the added 
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FIG. 4. Three-dimensional representation of the reduced den- 
sity of SPC water around a benzene molecule obtained by 
minimization of the functional. Blue to red indicate low to 
high densities up to p(r)/po — 3.5. The transparent grey sur- 
face that appears above the molecule represents the isosurface 
p(r)/p = 2.0. 
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FIG. 6. Site-oxygen radial distribution functions for a 
propane molecule in SPC water: MD results (red line) com- 
pared to the DFT results with and without the HS bridge term 
with R = 1.27 A (black and green lines, respectively). The top 
figure is a representation of the molecule and of the solvent 
reduced density isosurface corresponding to p(r)/po = 2.0. 




bridge term in this parameter range, see Figs [5] and |6j. 
Furthermore we show that using either the CS or PY 
version of the HS functional has a very small -although 
measurable- effect on the results. Retaining the CS ver- 
sion, we find that the MD value for methane in Table 
1 is reached when R ~ 1.27 A. Keeping that value, we 
find a close correlation to the MD results for the whole 
series of molecules. These encouraging findings are listed 
in Table 1 and further depicted in Fig. [8j Note that for 
each point in the figure, the computational effort to get 
the solvation free energy is orders of magnitude lower for 
3D-DFT than for MD. 



V. CONCLUSION 



FIG. 5. Site-oxygen radial distribution functions for a ben- 
zene molecule in SPC water: MD results (red line) compared 
to the DFT results with and without the HS bridge term with 
R — 1.27 A (black and green lines, respectively). 



This paper has presented the first three-dimensional 
implementation, to our knowledge, of the Kicrlik- 
Rosinberg fundamental measure theory for hard-sphere 
mixtures. Since the free-energy is written in terms of 
convolutions of the microscopic density with respect to 



Molecule 


MD 


DFT/HRF 


DFT/HRF+bridge 


methane 


10.96±0.46 


16.03 


10.43 


ethane 


10.75±0.50 


18.50 


10.33 


propane 


13.81±0.54 


24.86 


13.64 


butane 


14.69±0.54 


28.56 


14.71 


pentane 


15.43±0.59 


32.23 


15.87 


hexane 


16.40±0.63 


35.99 


17.02 


Ne 


11.21±0.46 


14.47 


11.67 


Ar 


8.661±0.46 


14.04 


9.51 


Kr 


8.242±0.46 


14.68 


9.20 



TABLE I. Solvation free energies of rare-gas atoms and alcane 
molecules in SPC water computed by DFT using the HRF 
approximation (eq. 18 with Tb = 0) or the HRF + hard- 
sphere bridge approximation (eqs [18] and eq. [19] with R = 
1.27A). They are compared to the MD values of RefP 3 for 
alcanes and RefP^for rare gases (with the corresponding error 
bars). All values are in kj/mol. 




1.20 



1.22 



1.24 1.26 
R[A] 



1.28 



1.30 



FIG. 7. Hydration free energy versus reference hard-sphere 
radius for methane, using either the Percus-Yevick (red circles 
and line) or Carnahan-Starling (black dashed lin e) versions of 
the FMT functional for the bridge term (eq. 19 1. 



discontinuous weight functions, the FMT density func- 
tionals are reputed to require very fine grids and very long 
recursive minimizations, at least for low dimensional ap- 
plications. In the three-dimensional case, we showed that 
using discrete 3D-FFT's which are well suited to describe 
the convolution of smooth functions with discontinuous 
(step-like or delta-like) distributions and using a more 
elaborated minimizer that the Picard iteration scheme 
that is usually prescribed^, the Kierlik-Rosinberg FMT 
functional can be efficiently minimized for realistic sys- 
tems with a grid resolution of only a few points per 
Angstrom and within at most a few tens of iterations. 
Such implementation constitutes a basis for tackling very 
diverse problems in physical chemistry involving fluids or 




FIG. 8. Comparison of the free energy of hydration for rare 
gases and alcanes (in trans conformation), calculated either 
by DFT with a hard-sphere bridge function of radius 1.27A 
or by MD simulations (Ref9^3 ancP^). From left to right: Kr, 
Ar, ethane, methane, propane, Ne, butane, pentane, hexane. 
The vertical bars correspond to the MD error bars indicated 
in Table 1. Units are kj/mol. 



solutions in the presence of atomistically-resolved inter- 
faces or confinement matrices. To describe realistic in- 
teractions, dispersive and Coulombic contributions can 
be easily introduced as m ean field perturbations to the 
hard-core functiona l 64 * 83 ^. We are working presently on 
several applications in that context. 

We have used here the implementation for another pur- 
pose: try to infer N-body corrections in the molecular 
density functional theory description of water, a system 
for which the restriction to two-body correlations (the 
so called homogeneous reference fluid approximation or 
HNC approximation in an inte gral e quation context) was 
found to present shortcoming d 14 * 34 ^. Hard-sphere correc- 
tions do seem to help for the special but important case 
of hydrophobic solvation. As stated in the introduc- 
tion, however, the present approach is meant to describe 
hydrophobicity at microscopic length scales but not to 
account properly for macroscopic phenomena such as 
dewettin gj 40 * 73 ! 74 !. Such limitation could be bypassed by 
adding coarse-grained contributions to the microscopic 
short-ranged functional. 

Furthermore the present effort has to be continued 
for the general situation of polar and charged solutes, 
for which electrostatic interactions and solvent angular 
dependence have to be included. In that case, despite 
positive results when M DFT is used as a post-treatment 
of exact densitied 36 * 37 *, we have some preliminary indica- 
tions, and there are premises in the literature toe) 6 * 15 * 89 -!, 
that life might not be so simple with self-consistent min- 
imization, and that angular-independent bridge correc- 
tions might not be sufficient. Mixing the type of hard- 
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sphere corrections studied here to the H-bonding three- 
body corrections introduced in RefP^ might be a way 
out. 
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Appendix A: A few useful formulas 

The first derivatives of the function &({n a }) in eq. [6] 
with respect to the weighted densities are required to 
compute the functional gradients and perform the mini- 
mization (see the Euler-Lagrange equation, eq. [9]). The 
second derivatives make it possible to compute the hard- 
sphere direct correlation functions, eq. |21| They are re- 
quired too if a N ewton-like minimization algorithm such 
as GMRE3M22lis usec j mst ead of L-BFGS (which is quite 
efficient but memory demanding). 

All those derivatives are listed below for the PY version 



of KR-FMT, eq. 16 



dn 


= - ln(l - na) 


d$ PY 


n 2 




dni 


1 - n 3 






'I" 1 + 


n\ 


dn 2 


1 - n 3 


87r(l-n 3 ) 2 


d$ PY 


n 


n\n 2 


dn 3 


1 - n 3 


(l-n 3 ) 2 1 


d 2 $ PY 


d 2 <s> PY 


1 


dnodn 3 


dn\dn 2 


1 - n 3 


d 2 $ PY 


n 2 




dn,idn 3 


(l-n 3 ) 2 






n 2 




dn\ 


47r(l — n 3 


) 2 


d 2 $ PY 


ni 


n 2 , 

+ 


dn 2 dn 3 


(l-n 3 ) 2 


47r(l - n 3 


d 2 $ PY 


n 


2n\n 2 




(l-n 3 ) 2 


+ (i-n 3 r 



12tt(1 - n 3 ) 3 



(Al) 



and for the CS version, eq. 16 

Q^CS 

ln(l-n 3 ) 

"2 



dn 

Q$CS 

dn\ 

q$CS 

dn 2 
g^CS 

dn 3 



l-n 3 



1 - n 3 12-7r(l - n 3 ) 2 n 3 
no — n 2 ! /(367m 3 ! ) nin 2 



12wnl 



ln(l - n 3 ) 



1 - n 3 



<9 2 $ 



cs 



187rn 3 (l - n 3 )3 
d 2$cs x 



(l-n 3 ) 2 36tt(1 - n 3 ) 2 n 2 
ln(l-n 3 ) 



187rn| 



dnodn 3 

Q2$CS 

dnidn 3 

Q2$CS 

dn 2 



dn±dn 2 

n 2 
(l-n 3 ) 2 



1 - n 3 



(A2) 



n 3 ) 



n 2 
6im 2 



ln(l - n 3 ) 



<9 2 $ 



cs 



dn 2 dn 3 



<9 2 $ 



cs 



dn 2 



4tt(1 - n 3 y 



67rn 3 (l 

— (n 3 (71-2(2 — 5n 3 + ri 3 ) — 127Trti(l — n 3 )n 2 ) 

+2n 2 (l - n 3 ) 3 ln(l - n 3 ))/(127r(l - n 3 fn 3 3 ) 

(n 3 (nl(6 - 21?i 3 + 26^ - 5n 3 ) + 727min 2 (l - n 3 )i 

+367rn (l - n 3 ) 2 n%\) + 6nf(l - n 3 f ln(l - n 3 )\ 
/(36tt(1 - n 3 ) 4 nj) 



All the second derivatives that are not written are equal 
to zero. There are thus 6 non-vanishing second deriva- 
tives to be considered instead of 21 in the Rosenfeld's 
vectorial version^. I n this respect also, the Kierlik- 
Rosinberg version of FMT appears much simpler to ma- 
nipulate and will be more efficient in Newton-like mini- 
mization schemes. 
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